function [J, grad] = lrCostFunction(theta, X, y, lambda)
    m = length(y);
    theta_reg = [0; theta(2:end)];
    J = -1 * sum(y .* log(sigmoid(X * theta)) + (1 - y) .* log(1 - sigmoid(X * theta))) / m + lambda / (2 * m) * theta_reg' * theta_reg;
    grad = X' * (sigmoid(X * theta) - y) / m + lambda / m * theta_reg;
end
